Дан график потребительских расходов в Австралии за период от 1989 до 1995 года по кварталам в единицах миллионов долларов. Требуется оценить дальнейшее поведение графика.
Загружаем данные и визуалилируем их:
Попробуем поделить данные на число дней в квартале с целью устранения шума
Ничего не изменилось, поэтому будем рассматривать исходные данные, без деления.
Построим STL-декомпозицию
Значения ошибок такие же больше по амплитуде, как и сезонность. Это вызвано тем, что дисперсия графика растет со временем, что хорошо видно на исходных данных. Чтобы исключить этот эффект, нужно применить преобразование Бокса-Кокса (и визуалилировать результат преобразования).
## [1] 0.1289189
Теперь дисперсия не увеличивается с изменением времени.
Сделаем ручной подбор модели ARIMA, для этого нужно сначала проверить стационарность данных
print(kpss.test(box_cox_transformed_series)$p.value)
## Warning in kpss.test(box_cox_transformed_series): p-value smaller than
## printed p-value
## [1] 0.01
Получаем, что ряд нестационарен. Применим сезонное дифференцирование:
Проверим снова на стационарность:
print(kpss.test(diff_series)$p.value)
## [1] 0.06700494
Конечно, p получилось больше 0.05, но все равно не сильно. На всякий случай лучше еще раз продифференцировать.
И снова проверяем стационарность.
## Warning in kpss.test(double_diff_series): p-value greater than printed p-
## value
## [1] 0.1
На этот раз p большое, можно считать, что ряд стационарен.
Построим ACF и PACF
Исходя из графиков, можно рассмотреть в качестве начального приближения модели ARIMA(1,1,1)(2,1,1)\(_4\). Рассмотрим несколько близких к ней моделей и посчитаем AICc для каждой:
| Модель | AICc |
|---|---|
| ARIMA(1,1,1)(2,1,1)\(_{4}\) | -491.3076563 |
| ARIMA(0,1,1)(2,1,1)\(_{4}\) | -490.0088317 |
| ARIMA(1,1,0)(2,1,1)\(_{4}\) | -489.9499495 |
| ARIMA(1,1,1)(2,1,0)\(_{4}\) | -479.9596413 |
| ARIMA(1,1,1)(1,1,1)\(_{4}\) | -491.3678759 |
| ARIMA(1,1,2)(2,1,1)\(_{4}\) | -489.13115 |
| ARIMA(1,1,1)(2,1,2)\(_{4}\) | -490.3129881 |
| ARIMA(1,1,1)(3,1,1)\(_{4}\) | -492.127228 |
| ARIMA(0,1,0)(2,1,1)\(_{4}\) | -491.4155876 |
| ARIMA(1,1,1)(1,1,0)\(_{4}\) | -459.3412136 |
| ARIMA(1,1,0)(1,1,1)\(_{4}\) | -490.0123476 |
| ARIMA(0,1,1)(1,1,1)\(_{4}\) | -490.0526418 |
| ARIMA(1,1,1)(0,1,1)\(_{4}\) | -493.3923725 |
| ARIMA(2,1,1)(0,1,1)\(_{4}\) | -491.2466169 |
| ARIMA(1,1,2)(0,1,1)\(_{4}\) | -491.2456675 |
| ARIMA(1,1,1)(0,1,2)\(_{4}\) | -491.4185055 |
Самая лучшая модель – ARIMA(1,1,1)(0,1,1)\(_4\)
Построим графики остатков
Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков
Тут все хорошо, теперь построим QQ-график и гистограмму остатков
Остатки похожи на нормальные. Но это, вместе со стационарностью и несмещенностью можно проверить с помощью критериев:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | не отвергается | 0.0590203 |
| Несмещённость | Уилкоксона | не отвергается | 0.849094 |
| Стационарность | KPSS | не отвергается | 0.1 |
Разобъем нашу выборку на обучающую и тестовую и посмотрим результат работы прогнозатора:
## ME RMSE MAE MPE MAPE
## Training set -17.17858 360.8368 263.0813 -0.04074299 0.7666702
## Test set -1006.41996 1144.9687 1006.4200 -1.76829683 1.7682968
## MASE ACF1 Theil's U
## Training set 0.2143337 -0.1205521 NA
## Test set 0.8199355 0.4816714 0.3549949
Хороший результат.
Обучаем автоматическую ARIMA и строим графики остатков
## Series: time_series
## ARIMA(1,0,0)(2,1,1)[4] with drift
## Box Cox transformation: lambda= 0.1289189
##
## Coefficients:
## ar1 sar1 sar2 sma1 drift
## 0.9127 -0.0594 -0.1584 -0.6406 0.0353
## s.e. 0.0402 0.1251 0.1084 0.1128 0.0027
##
## sigma^2 estimated as 0.001566: log likelihood=254.78
## AIC=-497.57 AICc=-496.94 BIC=-479.92
AICc слегка лучше, чем в ручной модели.
Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков
Тут все хорошо. Построим QQ-график и гистограмму остатков
Остатки похожи на нормальные. Проверим нормальность, несмещенность и стационарность критериями:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | не отвергается | 0.1601078 |
| Несмещённость | Уилкоксона | не отвергается | 0.8359361 |
| Стационарность | KPSS | не отвергается | 0.1 |
Все хорошо.
Разобъем выборку на обучающую и тестовую и посмотрим качество прогноза
## ME RMSE MAE MPE MAPE
## Training set -17.17858 360.8368 263.0813 -0.04074299 0.7666702
## Test set -1006.41996 1144.9687 1006.4200 -1.76829683 1.7682968
## MASE ACF1 Theil's U
## Training set 0.2143337 -0.1205521 NA
## Test set 0.8199355 0.4816714 0.3549949
Прогноз хороший.
Для ставнения двух полученных моделей сначала построим взаимный график их остатков
Теперь прогоним критерий Диболда-Мариано
##
## Diebold-Mariano Test
##
## data: resres.auto
## DM = 0.43677, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.663
## alternative hypothesis: two.sided
Критерий не говорит нам о том, что одна модель существенно лучше другой. Но, посмотря на AICc, на нормальность гистограмм остатков, видим, что все-таки автоматическая модель лучше.
Построим прогноз ETS
ets_model = ets(time_series, lambda=box_cox_lambda)
print(ets_model)
## ETS(A,A,A)
##
## Call:
## ets(y = time_series, lambda = box_cox_lambda)
##
## Box-Cox transformation: lambda= 0.1289
##
## Smoothing parameters:
## alpha = 0.7059
## beta = 1e-04
## gamma = 0.1704
##
## Initial states:
## l = 19.4323
## b = 0.0351
## s=-0.0182 -0.1311 0.1987 -0.0494
##
## sigma: 0.0394
##
## AIC AICc BIC
## -197.4138 -196.0705 -170.6855
Посмотрим на остатки:
Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков
Картина выгядит хуже, чем для аримы, но приемлима.
Построим QQ-график и гистограмму остатков
Остатки выглядят нормальными. Проверим нормальность, несмещенность и стационарность критериями:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | не отвергается | 0.119335 |
| Несмещённость | Уилкоксона | не отвергается | 0.808132 |
| Стационарность | KPSS | не отвергается | 0.1 |
Все хорошо.
Попробуем разбить временной ряд на обучающий и тестовый и посмотреть качество прогноза
## ME RMSE MAE MPE MAPE
## Training set -18.69817 362.1703 273.7808 -0.01541468 0.8115788
## Test set -1660.24570 1733.9249 1660.2457 -2.86729008 2.8672901
## MASE ACF1 Theil's U
## Training set 0.2230506 0.04478185 NA
## Test set 1.3526106 0.45730178 0.5329765
Хороший прогноз.
Сначала построим взаимный график остатков
Теперь попытаемся прогнать критерий Диболда-Мариано
##
## Diebold-Mariano Test
##
## data: res.autores.ets
## DM = -0.27691, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.7823
## alternative hypothesis: two.sided
Он не обнаруживает качественных отличий, однако, достигаемые уровни значимости критерия Льюнга-Бокса для остатков были лучше у аримы, и показатель AICc тоже. И, наконец, при разбиении на обучающую и тестовую временную серию, арима дала лучший результат.
Поэтому выбираем в качестве лучшей модели автоматическую ариму.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1995 Q4 63897.38 63104.15 64613.07 62455.42 65094.46
## 1996 Q1 67588.68 66495.47 68630.62 65899.40 69127.77
## 1996 Q2 62118.28 60966.45 63241.71 60216.40 63788.01
## 1996 Q3 64245.70 62952.68 65548.64 62171.56 66158.15
## 1996 Q4 65856.31 64275.27 67390.14 63380.12 68156.67
## 1997 Q1 69765.51 67939.59 71541.78 66894.58 72461.86
## 1997 Q2 64141.11 62346.52 65920.73 61301.19 66783.56
## 1997 Q3 66289.89 64341.89 68209.24 63323.91 69189.15
## 1997 Q4 68108.26 66013.48 70212.74 64816.16 71234.89
## 1998 Q1 72137.36 69799.97 74451.22 68504.67 75516.85
## 1998 Q2 66298.53 64055.68 68549.32 62917.96 69552.63
## 1998 Q3 68575.01 66174.31 70850.28 64954.50 72037.66
## 1998 Q4 70472.80 67902.15 72953.37 66558.60 74192.70
## 1999 Q1 74602.43 71829.34 77262.34 70258.35 78693.10
## 1999 Q2 68591.96 65900.90 71130.89 64555.09 72579.67
## 1999 Q3 70939.11 68075.44 73607.25 66666.48 75207.75
## 1999 Q4 72864.54 69871.68 75737.11 68320.76 77362.23
## 2000 Q1 77115.95 73868.06 80213.78 72110.21 81889.19
## 2000 Q2 70936.17 67849.12 73917.26 66180.91 75477.04
## 2000 Q3 73341.71 70082.74 76516.40 68406.04 78207.06
## 2000 Q4 75321.19 71856.37 78668.49 70092.75 80582.64
## 2001 Q1 79700.07 76022.15 83288.75 74064.07 85142.69
## 2001 Q2 73338.30 69874.32 76725.75 67981.31 78565.90
## 2001 Q3 75814.70 72138.40 79340.52 70319.03 81237.76
## 2001 Q4 77857.07 74003.56 81628.14 72066.86 83703.78
## 2002 Q1 82363.27 78253.63 86398.46 76246.74 88526.70
## 2002 Q2 75814.44 71896.76 79632.21 69944.86 81576.43
## 2002 Q3 78365.17 74316.33 82315.97 72228.45 84464.80
## 2002 Q4 80467.22 76231.64 84576.27 74146.77 86890.15
## 2003 Q1 85103.66 80615.51 89450.63 78331.31 91909.59